In [1]:
import requests
import json
import pprint
import numpy as np
import pandas as pd
from functools import reduce
from pandas.tools.plotting import scatter_matrix
from matplotlib import pyplot as plt
import itertools
import statsmodels.formula.api as smf
from sklearn.metrics import mean_squared_error
from math import sqrt
from IPython.core.display import display, HTML
import warnings
warnings.filterwarnings('ignore')
display(HTML("<style>.container { width:100% !important; }</style>"))
0 Introduction¶
What is happiness? This must be one of the most complicated questions to attack. In stead of giving any answer, this report examines data from newly released World Happniess Report 2019, integrate it with other data from World Bank API and explores various combinations of indicators. With such data being explored, linear regression models are fit and the results are studied. At last, a mutiple-variable linear model is presented and interpreted, with which happniess at least can be predicated, even though we still have no idea what happniess is.
1 Data wrangling¶
1.a Static Data¶
In [2]:
wh = pd.read_excel('WHRData.xls', index_col=0)
wh = wh.reset_index()
wh = wh.rename(columns={'Country name': 'country', 'Year': 'year'})
wh.head()
Out[2]:
1.b Dynamic data (using world bank api)¶
In [3]:
class WorldBankDataFetcher:
"""
"""
def __init__(self, indicators):
self.indicators = indicators
self.jsons = self.dataframes = []
self.merged_data = None
def fetch_json_by(self, indicator):
url = "http://api.worldbank.org/v2/countries/all/indicators/%s&format=json&per_page=9000" % (indicator + '?date=2000:2018')
r = requests.get(url)
return json.loads(r.content)
def json_to_df(self, j):
indicator_name = j[1][0]['indicator']['value']
df = pd.DataFrame.from_dict(j[1])
df['country'] = df[['country']].applymap(lambda x : x['value'])
df = df[['country', 'countryiso3code', 'date', 'value']]
df.columns = ['country', 'code', 'year', indicator_name]
return df
def fetch(self):
print('Note: this operation make take a while depending on the network status...\n...')
for i in indicators:
self.jsons.append(self.fetch_json_by(i))
progress = str(indicators.index(i)+1) + '/' + str(len(indicators))
print(' '.join([progress, 'Indicator', i]))
print('Creating data frames...')
self.dataframes = [self.json_to_df(j) for j in self.jsons]
print('Merging data frames...')
self.merged_data = reduce(lambda x, y: x.merge(y, how='outer', on=['country', 'code', 'year']), self.dataframes)
print('All done!')
Indicators choosing process¶
- The featured indicator list can be found here: https://data.worldbank.org/indicator
- By clicking any of the indicator, a basic chart for this indicator can be found
- Another priciple I followed is choose indicators that come from different sectors
- The chosen indicators was listed as below
In [4]:
indicators = [
'AG.LND.FRST.K2', # Forest area (sq. km)
'AG.LND.AGRI.ZS', # Agricultural land (% of land area)
'EG.ELC.ACCS.ZS', # Access to electricity (% of population)
'EN.ATM.CO2E.PC', # CO2 emissions (metric tons per capita)
'MS.MIL.XPND.GD.ZS', # Military expenditure (% of GDP)
'SL.UEM.TOTL.ZS', # Unemployment, total (% of total labor force) (modeled ILO estimate)
'SE.PRM.UNER.FE', # Children out of school, primary, female
'SP.POP.TOTL', # Population, total
'IS.AIR.DPRT', # Air transport, registered carrier departures worldwide
'IC.TAX.TOTL.CP.ZS', # Total tax rate (% of commercial profits)
'NY.GDP.MKTP.KD.ZG', # GDP growth (annual %)
'SE.ADT.LITR.ZS' # Literacy rate, adult total (% of people ages 15 and above)
]
In [5]:
fetcher = WorldBankDataFetcher(indicators)
fetcher.fetch()
In [6]:
wb = fetcher.merged_data
wb['year'] = pd.to_numeric(wb['year'])
wb = wb[wb.code!='']
wb = wb.drop('code', axis=1)
wb.head()
Out[6]:
1.c Integration¶
- Construct a contry name dictionary by looking at the set differnence between the two data frames
- Use the contry name dictionary to make the country name consisteny between the two data frames
- Merge the two data frames into a single one
In [7]:
set(wh.country).difference(set(wb.country))
Out[7]:
In [8]:
set(wb.country).difference(set(wh.country))
Out[8]:
In [9]:
country_names = {
'Congo (Brazzaville)' : 'Congo, Rep.',
'Congo (Kinshasa)' : 'Congo, Dem. Rep.',
'Egypt' : 'Egypt, Arab Rep.',
'Gambia': 'Gambia, The',
'Hong Kong S.A.R. of China' : 'Hong Kong SAR, China',
'Iran' : 'Iran, Islamic Rep.',
'Ivory Coast' : "Cote d'Ivoire",
'Kyrgyzstan' : 'Kyrgyz Republic',
'Laos' : 'Lao PDR',
'Macedonia': 'North Macedonia',
'Palestinian Territories' : 'West Bank and Gaza',
'Russia' : 'Russian Federation',
'Slovakia' : 'Slovak Republic',
'South Korea' : 'Korea, Rep.',
'Swaziland': 'Eswatini',
'Syria' : 'Syrian Arab Republic',
'Venezuela' : 'Venezuela, RB',
'Yemen' : 'Yemen, Rep.'
}
country_names = {v:k for k, v in country_names.items()}
country_names
Out[9]:
In [10]:
wb = wb.applymap(lambda x: country_names[x] if x in country_names else x)
In [11]:
set(wh.country).difference(set(wb.country))
Out[11]:
Pre-integration check¶
In [12]:
len(wh)
Out[12]:
Integration¶
In [13]:
integrated = wh.merge(wb, how='left', on=['country', 'year'])
Post-integration check¶
In [14]:
len(integrated)
Out[14]:
In [15]:
integrated.head()
Out[15]:
- Since the row count of the integrated data frame is equae to the World Happniess Report data frame, hence no data point was lost duing this mergeing.
2 EDA¶
2.a EDA – data quality inspection¶
- Check null value count for each column
- Check and drop, if necessary, duplicate rows
- Drop columns with less than 1500 non-null values
- Drop rows with any null values
- Check target variable distribution and transform it if necessary
In [16]:
integrated.isnull().sum().sum()
Out[16]:
In [18]:
integrated.isnull().sum()
Out[18]:
In [19]:
integrated.duplicated().any()
Out[19]:
In [20]:
dropped = integrated.dropna(axis='columns', thresh=1500)
dropped = dropped.dropna(axis='index', thresh=len(dropped.columns))
dropped.isnull().sum().sum()
Out[20]:
In [21]:
dropped.isnull().sum()
Out[21]:
In [22]:
dropped = dropped.reset_index(drop=True)
dropped = dropped.loc[:, dropped.columns != 'country']
dropped.head()
Out[22]:
In [23]:
data = dropped.copy()
data['year'] = data['year'] - 2000
_ = data[['Life Ladder']].hist(bins=50, figsize=(5,5), alpha=0.8)
plt.show()
- As one of the assumptions that Linear Regresion model need meet, target variable should be normal distribution.
- As shown above, the distribution of the indicator Life Ladder basicly is a normal distribution. Thus, it is not necessary to transform it.
2.b EDA – the search for good predictors (12 marks)¶
- Rename columns such that the data set is better for visualization
- Draw histogram plot for each feature and choose features that comform to normal distribution
- Draw scatter matrix and correlation matrix to choose features that are correlated to happniess and not quite correlated between them.
In [24]:
name_table = {
'Life Ladder': 'happiness',
'Log GDP per capita': 'gdp_pc',
'Social support': 'support',
'Population, total': 'pop',
'Total tax and contribution rate (% of profit)': 'tax',
'GDP growth (annual %)': 'gdp_g',
'Healthy life expectancy at birth': 'life_expected',
'Freedom to make life choices': 'choice',
'Generosity': 'generosity',
'Perceptions of corruption': 'corruption',
'Positive affect': 'positive',
'Negative affect': 'negative',
'Confidence in national government': 'gov_confidence',
'Democratic Quality': 'democoratic',
'Delivery Quality': 'delivery',
'Standard deviation of ladder by country-year': 'sd_country_year',
'Standard deviation/Mean of ladder by country-year': 'sd_m_country_year',
'GINI index (World Bank estimate), average 2000-16': 'gini',
'Unemployment, total (% of total labor force) (modeled ILO estimate)': 'unemployment'
}
data = data.rename(columns=name_table)
data.head()
Out[24]:
In [25]:
f, ax = plt.subplots(2, 10, figsize=(30, 10))
for i in range(0, len(data.columns)):
c = 10
f = data.columns[i]
sub = ax[i//c, i%c]
sub.hist(data[f], bins=30, alpha=0.8)
sub.set_title(f)
plt.show()
- Feartures year, pop and tax are dropped since they are extremely skewed
In [26]:
_ = scatter_matrix(data[data.columns.difference(['pop', 'tax', 'year'])], alpha=0.6, figsize=(32, 32))
plt.show()
In [27]:
data[data.columns.difference(['pop', 'tax', 'year'])].corr()
Out[27]:
- As shown above, features gov_confidence and unemployment have low correlations to happniess. Thus, the two feature should not be chosen.
- All three features , life_expected, delivery , democoratic and gdp_pc indicate high correlations to happniess. However, the correlations among them are also very high. To preven multicollinearity, only one should be chosen
- sd_country_year and sd_m_country_year are hard to intrepret, thus they will not be considered as well.
- Based on the above the chosen features are listed below:
In [28]:
chosen_features = [
'gdp_pc',
'choice',
'generosity',
'corruption',
'positive',
'negative',
'gini',
'gdp_g',
'support'
]
In [29]:
all_combinations = []
for length in range(1, len(chosen_features)+1):
for subset in itertools.combinations(chosen_features, length):
all_combinations.append(subset)
all_combinations
Out[29]:
3.2 Fit model for each of the combinations and record values such as F test, AI, RMSE, etc¶
In [30]:
all_combinations_df = pd.DataFrame()
for c in all_combinations:
all_combinations_df = all_combinations_df.append({f: 1 if f in c else 0 for f in chosen_features}, ignore_index=True)
all_combinations_df['f test'] = all_combinations_df['r^2'] = all_combinations_df['aic'] = all_combinations_df['f_pvalue'] = all_combinations_df['rmse'] = all_combinations_df['no of p>0.05'] = all_combinations_df['df'] = 0
all_combinations_df['formula'] = ''
all_combinations_df.head()
Out[30]:
In [31]:
for i, row in all_combinations_df.iterrows():
features = [f for f in chosen_features if row[f] == 1]
formula = 'happiness ~ 1 + ' + ' + '.join(features)
mod = smf.ols(formula=formula, data=data).fit()
all_combinations_df.iloc[i, all_combinations_df.columns.get_loc("formula")] = formula
all_combinations_df.iloc[i, all_combinations_df.columns.get_loc("aic")] = mod.aic
all_combinations_df.iloc[i, all_combinations_df.columns.get_loc("f test")] = mod.fvalue
all_combinations_df.iloc[i, all_combinations_df.columns.get_loc("r^2")] = mod.rsquared
all_combinations_df.iloc[i, all_combinations_df.columns.get_loc("f_pvalue")] = mod.f_pvalue
all_combinations_df.iloc[i, all_combinations_df.columns.get_loc("rmse")] = sqrt(mean_squared_error(data.happiness, mod.fittedvalues))
all_combinations_df.iloc[i, all_combinations_df.columns.get_loc("no of p>0.05")] = len([p for p in mod.pvalues if p > 0.05])
all_combinations_df.iloc[i, all_combinations_df.columns.get_loc("df")] = mod.df_model
all_combinations_df.head()
Out[31]:
3.3 Filtter out models with parameters that have p value greater than 0.05¶
In [32]:
all_significant = all_combinations_df[all_combinations_df['no of p>0.05'] == 0]
all_significant
Out[32]:
3.4 Narrow down to models with high R-Squared value¶
In [33]:
all_significant_high_rsq = all_significant[all_significant['r^2']>0.76]
all_significant_high_rsq
Out[33]:
3.5 Narrow down to models with low RMSE value¶
In [34]:
all_significant_high_rsq_low_rmse = all_significant_high_rsq[all_significant_high_rsq.rmse < 0.56]
all_significant_high_rsq_low_rmse
Out[34]:
3.6 Final model¶
* According to okam razors, the model with index=384 is chosen, as shown below.In [35]:
final_model = smf.ols(formula='happiness ~ 1 + choice + corruption + gdp_pc + generosity + positive + support', data=data).fit()
3.7 Model interpretation¶
In [36]:
x = final_model.fittedvalues
y = final_model.resid
plt.figure(figsize=(12, 12))
plt.scatter(x, y, c="g", alpha=0.4, s=18)
plt.xlabel("Fitted Values")
plt.ylabel("Residuals")
plt.title('Residuals vs Fitted Values')
plt.plot(np.unique(x), np.poly1d(np.polyfit(x, y, 1))(np.unique(x)), c='b')
plt.show()
In [37]:
final_model.summary()
Out[37]:
Residuals:¶
- As shown above no obvious pattern or structure can be found. Hence, the homoscedasticity assumption is met.
- Also the residuals are randomly scatted around the fitted line. Hence the linearity assumption is met.
Equaltion, represented as vector product: [choice, corruption, gdp_pc, generosity, positive, support] * [0.55, -0.87, 0.49, 0.51, 1.90, 2.13].transpose()¶
- For every 1 unit increase in choce, happiness rises by 0.55
- For every 1 unit decrease in corruption, happiness rises by 0.87
- For every 1 unit increase in gdp_pc, happiness rises by 0.49
- For every 1 unit increase in generosity, happiness rises by 0.51
- For every 1 unit increase in positive, happiness rises by 1.90
- For every 1 unit increase in support, happiness rises by 2.13
Reliability¶
- R-squared = 0.77 makes this model a good fit.
- With all assuptions being met, this model is worth recommending this model to a client
4 Conclusion¶
Happniess can be predicated with various models. As disccussed above, the best combination of indicators for this purpose is Positive affect, Social support, GDP per capita, Freedom to make life choices, Generosity and Corruption. The linear regression model built with this combination show a good fit with a R-Squared=0.77 and a RMSE=0.56.